
if !d.name eq 'PS' then begin
    device,xsize=18,ysize=12,yoffset=3
    !p.thick=2 & !x.thick=1 & !y.thick=1
end

@eu_setup
@rhoprof
restore,'urms.sav'
nn=where(r ge 0.72)
urms=mean(surms[nn])
; density heigth scale
hd=-1./xder_curl(alog(rh),r)
tau=Hd/urms

etat=10^4*rr*tau*surms^2/3.

rsinth1=fltarr(m,l)
dwdphi=fltarr(n,m,l,ns)
dvdphi=fltarr(n,m,l,ns)
dsinudth=fltarr(n,m,l,ns)
dwdth=fltarr(n,m,l,ns)
drvdr=fltarr(n,m,l,ns)
drudr=fltarr(n,m,l,ns)

wr=fltarr(n,m,l,ns)
wth=fltarr(n,m,l,ns)
wphi=fltarr(n,m,l,ns)

;TOTAL HELICITY COMPONENT

for it=0,ns-1 do begin
print,'Computing total helicity',it
 for k=0, l-1 do begin
  for i=0,n-1 do begin
   dsinudth[i,*,k,it]=xder_curl(reform(smooth(sin(theta)*u[i,*,k,it],3)),theta)       
   dwdth[i,*,k,it]=xder_curl(reform(smooth(w[i,*,k,it],3)),theta)
  endfor
 endfor

 for j=0,m-1 do begin
  for i=0,n-1 do begin
   drvdr[i,j,*,it]=xder_curl(smooth(r[*]*v[i,j,*,it],3),r)   
   drudr[i,j,*,it]=xder_curl(smooth(r[*]*u[i,j,*,it],3),r)
  endfor
 endfor
;

 for k=0, l-1 do begin
  for j=1,m-1 do begin
   rsinth1[j,k]=1./(r[k]*sin(theta[j]))
  endfor
 endfor
 rsinth1[0,*]=0.
 rsinth1[m-1,*]=0.


 for k=0, l-1 do begin
  for j=1,m-1 do begin
   for i=0,n-1 do begin
     wr[i,j,k,it]=rsinth1[j,k]*dsinudth[i,j,k,it]
     wth[i,j,k,it]=rsinth1[j,k]*dwdphi[i,j,k,it]-drudr[i,j,k,it]/r[k]
     wphi[i,j,k,it]=drvdr[i,j,k,it]/r[k]-dwdth[i,j,k,it]/r[k]
   endfor
  endfor
 endfor

endfor

Ht=total((wr*w + wth*v + wphi*u),1)/double(n)
Htm=total(Ht,3)/double(ns)

; CURL OF THE MEAN
print,'Computing mean helicity',ns
uum=total(u,1)/double(n)
vvm=total(v,1)/double(n)
wwm=total(w,1)/double(n)

um=total(uum,3)/double(ns)
vm=total(vvm,3)/double(ns)
wm=total(wwm,3)/double(ns)

dwdphim=fltarr(m,l)
dvdphim=fltarr(m,l)
dsinudthm=fltarr(m,l)
dwdthm=fltarr(m,l)
drvdrm=fltarr(m,l)
drudrm=fltarr(m,l)

wrm=fltarr(m,l)
wthm=fltarr(m,l)
wphim=fltarr(m,l)

 for k=0, l-1 do begin
   dsinudthm[*,k]=xder_curl(reform(smooth(sin(theta)*um[*,k],2)),theta)       
   dwdthm[*,k]=xder_curl(reform(smooth(wm[*,k],2)),theta)
 endfor
 for j=1,m-1 do begin
   drvdrm[j,*]=xder_curl(smooth(r[*]*vm[j,*],2),r)   
   drudrm[j,*]=xder_curl(smooth(r[*]*um[j,*],2),r)
 endfor

 for k=0, l-1 do begin
  for j=1,m-1 do begin
     wrm[j,k]=rsinth1[j,k]*dsinudthm[j,k]
     wthm[j,k]=rsinth1[j,k]*dwdphim[j,k]-drudrm[j,k]/r[k]
     wphim[j,k]=drvdrm[j,k]/r[k]-dwdthm[j,k]/r[k]
  endfor
 endfor
hm=wrm*wm + wthm*vm + wphim*um
khel=htm-hm
 for k=0, l-1 do begin
  for j=1,m-1 do begin
    khel[j,k]=-tau[k]*khel[j,k]/3.
 endfor
 endfor
save,filename='alpha_eta.sav',khel,etat
!p.multi=[0,2,1] & !p.charsize=1.3 
contour,transpose(smooth(khel,3)),x,y,/fi,nl=64,/iso,title='!4a!8=-!4s!8<!4x!8`.u`>/3!6'
plot,r,smooth(etat,3),/ylog,xtitle='!8r!6',ytitle='!4g!8!Dt!N!6(cm!U!82!N!6s!U!8-1!N!6)!6',$
       title='!4g!8!Dt!N=!4s!8u!D!6rms!8!U2!N/3!6',thick=3
end
